UPVacina

Bibliotecas

install.packages("qqplotr")

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into ‘C:/Users/Wasim/AppData/Local/R/win-library/4.3’
(as ‘lib’ is unspecified)
also installing the dependencies ‘DEoptimR’, ‘pracma’, ‘twosamples’, ‘pbmcapply’, ‘robustbase’, ‘opdisDownsampling’, ‘qqconf’

tentando a URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/DEoptimR_1.1-3.zip'
Content type 'application/zip' length 74909 bytes (73 KB)
downloaded 73 KB

tentando a URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/pracma_2.4.4.zip'
Content type 'application/zip' length 1733855 bytes (1.7 MB)
downloaded 1.7 MB

tentando a URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/twosamples_2.0.1.zip'
Content type 'application/zip' length 278308 bytes (271 KB)
downloaded 271 KB

tentando a URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/pbmcapply_1.5.1.zip'
Content type 'application/zip' length 47067 bytes (45 KB)
downloaded 45 KB

tentando a URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/robustbase_0.99-3.zip'
Content type 'application/zip' length 3081839 bytes (2.9 MB)
downloaded 2.9 MB

tentando a URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/opdisDownsampling_1.0.1.zip'
Content type 'application/zip' length 3473242 bytes (3.3 MB)
downloaded 3.3 MB

tentando a URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/qqconf_1.3.2.zip'
Content type 'application/zip' length 1420046 bytes (1.4 MB)
downloaded 1.4 MB

tentando a URL 'https://cran.rstudio.com/bin/windows/contrib/4.3/qqplotr_0.0.6.zip'
Content type 'application/zip' length 1111731 bytes (1.1 MB)
downloaded 1.1 MB

package ‘DEoptimR’ successfully unpacked and MD5 sums checked
package ‘pracma’ successfully unpacked and MD5 sums checked
package ‘twosamples’ successfully unpacked and MD5 sums checked
package ‘pbmcapply’ successfully unpacked and MD5 sums checked
package ‘robustbase’ successfully unpacked and MD5 sums checked
package ‘opdisDownsampling’ successfully unpacked and MD5 sums checked
package ‘qqconf’ successfully unpacked and MD5 sums checked
package ‘qqplotr’ successfully unpacked and MD5 sums checked

The downloaded binary packages are in
    C:\Users\Wasim\AppData\Local\Temp\RtmpQpORmZ\downloaded_packages

install.packages("qqplotr")

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/
Installing package into ‘C:/Users/Wasim/AppData/Local/R/win-library/4.3’
(as ‘lib’ is unspecified)

Unindo tabelas

Dados de estados

lista_UF_cobertura = list.files() 

# Lista todos os arquivos CSV no diretório
arquivos_xlsx <- list.files(here(), pattern = "\\.xlsx$", full.names = F)

#Criar funcao
read_xl_sipni <- function(arquivo) {
  dados <- read_excel(arquivo, skip = 4, col_names = FALSE) %>%
    setNames(.[1, ]) %>%
    slice(-1:-4) %>%
    mutate(across(-1, as.numeric)) %>%  
    pivot_longer(cols = -c(1), 
                 names_to = "imuno", 
                 values_to = "cobertura") %>% 
    mutate(arquivo = basename(arquivo)) %>% 
    janitor::clean_names() %>% 
    mutate(cobertura = round(cobertura, 2),
           ano = as.numeric(gsub(".*_(\\d{4})\\.xlsx", "\\1", arquivo)),
           uf = gsub("\\d", "", .[[1]]))  %>% 
    select(-unidade_da_federacao) %>% 
    filter(!is.na(uf),
           !grepl("Gerado", uf)) %>% 
    select(uf, ano, imuno, cobertura) %>% 
    group_by(uf, ano, imuno)
  
  return(dados)
}


# Leia, converta em tabela longa e combine os arquivos CSV em um único DataFrame
sipni_cobertura_uf_1994_2023 <- lapply(arquivos_xlsx, read_xl_sipni) %>%
  bind_rows()

saveRDS(sipni_cobertura_uf_1994_2023, file = "sipni_cobertura_uf_1994_2023.rds")

Dados de municípios

#Listar tabelas
lista_MU_cobertura = list.files() 
arquivos_xlsx <- list.files(here(), pattern = "\\.xlsx$", full.names = F)

#Criar função
read_xl_sipni_municipio <- function(arquivo) {
  dados <- read_excel(arquivo, skip = 4, col_names = FALSE) %>%
    setNames(.[1, ]) %>%
    slice(-1:-4) %>%
    mutate(across(-1, as.numeric)) %>%  
    pivot_longer(cols = -c(1), 
                 names_to = "imuno", 
                 values_to = "cobertura") %>% 
    mutate(arquivo = basename(arquivo)) %>% 
    janitor::clean_names() %>% 
    mutate(cobertura = round(cobertura, 2),
           ano = as.numeric(gsub(".*_(\\d{4})\\.xlsx", "\\1", arquivo))) %>% 
    filter(!is.na(municipio),
           !grepl("Gerado", municipio)) %>% 
    select(municipio, ano, imuno, cobertura) %>% 
    group_by(municipio, ano, imuno)
  
  return(dados)
}

#Unir dados
sipni_cobertura_municipios_1994_2023 <- lapply(arquivos_xlsx, read_xl_sipni_municipio) %>%
  bind_rows()

#Salvar
saveRDS(sipni_cobertura_municipios_1994_2023, 
          file = "sipni_cobertura_municipios_1994_2023.rds")

Processamento de dados

#Estados
sipni_cobertura_uf_1994_2023_2 = sipni_cobertura_uf_1994_2023 %>% 
  mutate(mun_uf = "UF", 
         nome = uf,
         nome =  toupper(nome),
         nome = stri_trans_general(nome, "Latin-ASCII")) %>% 
  select(nome, uf, ano, imuno, cobertura, mun_uf)

saveRDS(sipni_cobertura_uf_1994_2023_2, file = "sipni_cobertura_uf_1994_2023_2.rds")

#Municípios
sipni_cobertura_municipios_1994_2023_2 = sipni_cobertura_municipios_1994_2023 %>% 
 mutate(codigo = as.character(str_extract(municipio, "\\d+")),# Extrair números
        nome = str_remove(municipio, "\\d+ "),# Extrair texto
        mun_uf = "Município" ) %>% 
  select(!"...1":municipio)

glimpse(sipni_cobertura_municipios_1994_2023_2)

# #Unir
# sipni_all = bind_rows(sipni_cobertura_uf_1994_2023_2, sipni_cobertura_municipios_1994_2023_2)
# write.csv(sipni_all, file = "sipni_uf_mun_1994_2023.csv")

#Anotar municipios e estados

#Anotações de cidades. Fonte: IBGE. https://www.ibge.gov.br/explica/codigos-dos-municipios.php
dtb_municipios_cod = RELATORIO_DTB_BRASIL_MUNICIPIO %>% 
  clean_names() %>% 
  select(codigo_uf = uf, uf = nome_uf, codigo = codigo_municipio_completo, nome_municipio) %>% 
  mutate(nome_municipio_original = nome_municipio,
         nome =  toupper(nome_municipio_original),
         nome = stri_trans_general(nome, "Latin-ASCII")) %>% 
  select(-nome_municipio)

municipios = sipni_cobertura_municipios_1994_2023_2 %>% 
  mutate(nome = gsub("\\\\", "", nome)) %>% 
  left_join(dtb_municipios_cod, by = "nome") %>% 
  select(nome, nome_municipio_original, uf, codigo_municipio = codigo.y, ano, imuno, cobertura, codigo_uf, mun_uf)

#Salvar
saveRDS(municipios, file = "sipni_cobertura_municipios_1994_2023_2.rds")
  

#Análise de dados Os dados das populações foram obtidos dos anos 2000 a 2022, pois os anos 2007 e 2023 não estavam disponíveis. Definir variáveis - https://www.ibge.gov.br/estatisticas/downloads-estatisticas.html - Censos - Perfil Estados - Perfil Municipios - Economicos - Indicadores sociais - Educacação e qualificação profissional - Economia de saúde - Acesso à internet - Educação - Saúde

install.packages("basedosdados")
library("basedosdados")

# Para carregar o dado direto no R
query <- bdplyr("br_ibge_censo_demografico.microdados_domicilio_1970")
df <- bd_collect(query)


install.packages("bdplyr")
query <- bdplyr("br_ms_atencao_basica.municipio")
df <- bd_collect(query)

# População ---- 

# Obter população de municipios.
anos = 2000:2023
resultados = list()
for (ano in anos) {
  tryCatch({
    # Chamar a função populacao_municipios para o ano atual e armazenar o resultado na lista
    df_ano <- populacao_municipios(ano)
    df_ano <- df_ano %>%
      mutate_all(as.character)  # Convertendo todas as colunas para character
    resultados[[as.character(ano)]] <- df_ano
  }, error = function(e) {
    # Tratar o erro (por exemplo, imprimir uma mensagem)
    print(paste("Erro para o ano", ano, ":", conditionMessage(e)))
  })
}

# Combinar todos os dataframes em um único dataframe
populacao_municipios_2000_2022 <- bind_rows(resultados, .id = "ano") %>% 
  select(uf_abrev = uf, 
         nome_municipio_original = nome_munic,
         ano, codigo_uf, 
         populacao, 
         codigo_municipio = cod_municipio) %>% 
  mutate(ano = as.numeric(ano),
         populacao = as.numeric(populacao))
  
print(populacao_municipios_2000_2022)


# GDP ----

anos = 1999:2020
pib_municipios(ano = 2003)

resultados = list()
for (ano in anos) {
  tryCatch({
    df_ano <- pib_municipios(ano = ano, dir=".")
    resultados[[ano]] <- df_ano
  }, error = function(e) {
    print(paste("Erro para o ano", ano, ":", conditionMessage(e)))
  })
}

# Combinar todos os dataframes em um único dataframe
populacao_municipios_2000_2022 <- bind_rows(resultados, .id = "ano") %>% 
  select(uf_abrev = uf, 
         nome_municipio_original = nome_munic,
         ano, codigo_uf, 
         populacao, 
         codigo_municipio = cod_municipio) %>% 
  mutate(ano = as.numeric(ano),
         populacao = as.numeric(populacao))
  
print(populacao_municipios_2000_2022)

#Unir dados

municipios_2 = municipios %>% 
 left_join(populacao_municipios_2000_2022 %>% 
             select(codigo_municipio, ano, populacao),
            by = c("codigo_municipio", "ano"))

#IEPS dados

#Unir dados
IEPS_Brasil_2010_2022_Todos <- read_excel("dados_brutos/IEPS/IEPS_Brasil_2010_2022_Todos.xlsx") %>% 
  mutate(nivel = "Brasil")
IEPS_Estados_2010_2022_Todos <- read_excel("dados_brutos/IEPS/IEPS_Estados_2010_2022_Todos.xlsx") %>% 
  mutate(nivel = "Estado")
IEPS_MacroRegiao_2010_2022_Todos <- read_excel("dados_brutos/IEPS/IEPS_MacroRegiao_2010_2022_Todos.xlsx") %>% 
  mutate(nivel = "Macro Região")
IEPS_Municipios_2010_2022_Todos <- read_excel("dados_brutos/IEPS/IEPS_Municipios_2010_2022_Todos.xlsx") %>% 
  mutate(nivel = "Município")
IEPS_Regiao_2010_2022_Todos <- read_excel("dados_brutos/IEPS/IEPS_Regiao_2010_2022_Todos.xlsx") %>% 
  mutate(nivel = "Região")

IEPS_Completo_2010_2022_Todos = bind_rows(IEPS_Brasil_2010_2022_Todos,
                                          IEPS_Estados_2010_2022_Todos,
                                          IEPS_MacroRegiao_2010_2022_Todos,
                                          IEPS_Regiao_2010_2022_Todos,
                                          IEPS_Municipios_2010_2022_Todos) 

#Padronizar tabelas

IEPS_Brasil_2010_2022_Todos <- IEPS_Brasil_2010_2022_Todos %>% 
  mutate(nivel = "Brasil") %>% 
  mutate(across(c(cob_ab:pct_pop_masc), ~ str_replace(., ",", "."))) %>% 
  mutate(across(c(cob_ab:pct_pop_masc), ~ as.numeric(.))) %>% 
  mutate(ano = lubridate::ymd(ano, truncated = 2L))

IEPS_Estados_2010_2022_Todos <- IEPS_Estados_2010_2022_Todos %>% 
  mutate(nivel = "Estado") %>% 
  mutate(across(c(cob_ab:pct_pop_masc), ~ str_replace(., ",", "."))) %>% 
  mutate(across(c(cob_ab:pct_pop_masc), ~ as.numeric(.))) %>% 
  mutate(ano = lubridate::ymd(ano, truncated = 2L))

IEPS_MacroRegiao_2010_2022_Todos <- IEPS_MacroRegiao_2010_2022_Todos %>% 
  mutate(nivel = "Macro Região") %>% 
  mutate(across(c(cob_ab:pct_pop_masc), ~ str_replace(., ",", "."))) %>% 
  mutate(across(c(cob_ab:pct_pop_masc), ~ as.numeric(.))) %>% 
  mutate(ano = lubridate::ymd(ano, truncated = 2L))

IEPS_Regiao_2010_2022_Todos <- IEPS_Regiao_2010_2022_Todos %>% 
  mutate(nivel = "Região") %>% 
  mutate(across(c(cob_ab:pct_pop_masc), ~ str_replace(., ",", "."))) %>% 
  mutate(across(c(cob_ab:pct_pop_masc), ~ as.numeric(.))) %>% 
  mutate(ano = lubridate::ymd(ano, truncated = 2L))

IEPS_Municipios_2010_2022_Todos <- IEPS_Municipios_2010_2022_Todos %>% 
  mutate(nivel = "Município",
         porte_municipio = case_when(pop <= 20000 ~ "Município de Pequeno Porte I",
                                     pop <= 50000 ~ "Município de Pequeno Porte II",
                                     pop <= 100000 ~ "Município de Médio Porte",
                                     pop <= 900000 ~ "Município de Grande Porte",
                                     TRUE ~ "Metrópole")) %>% 
  mutate(across(c(cob_ab:pct_pop_masc), ~ str_replace(., ",", "."))) %>% 
  mutate(across(c(cob_ab:pct_pop_masc), ~ as.numeric(.))) %>% 
  mutate(ano = lubridate::ymd(ano, truncated = 2L))

IEPS_Completo_2010_2022_Todos = IEPS_Completo_2010_2022_Todos %>% 
  mutate(across(c(cob_ab:pct_pop_masc), ~ str_replace(., ",", "."))) %>% 
  mutate(across(c(cob_ab:pct_pop_masc), ~ as.numeric(.))) %>% 
  mutate(ano = lubridate::ymd(ano, truncated = 2L))

IEPS_Completo_2010_2022_Todos %>% 
  saveRDS(file = here("dados_processados", "IEPS_Completo_2010_2022_Todos.rds"))

Seleção de variáveis

#Análise exploratória

#Definir tema
ggthemr("fresh", spacing = 1)

Brasil

Estados

São Paulo

#Importar dados
IEPS_Completo_2010_2022_selecionados_SP <- readRDS(here("dados_processados", "IEPS_Completo_2010_2022_selecionados_SP.rds")) 
IEPS_var_selecionadas = read_excel("dados_processados/IEPS_codebook_selecionados.xlsx") %>% 
  clean_names() %>% 
  filter(selecao == "Sim") %>% 
  mutate(bloco = if_else(str_detect(variavel, "cob_vac"), "Cobertura vacinal", bloco))

#Manipulação
IEPS_Completo_filtrado_SP = IEPS_Completo_2010_2022_selecionados_SP %>% 
  filter(nivel == "Estado") %>% 
  pivot_longer(cols = cob_ab:pct_pop_masc,
               names_to = "variavel",
               values_to = "valor") %>% 
  inner_join(IEPS_var_selecionadas, by = "variavel") %>% 
  mutate(ano = lubridate::ymd(ano, truncated = 2L))


IEPS_Completo_Regiao_SP = IEPS_Completo_2010_2022_selecionados_SP %>% 
  filter(nivel == "Região") %>% 
  pivot_longer(cols = cob_ab:pct_pop_masc,
               names_to = "variavel",
               values_to = "valor") %>% 
  inner_join(IEPS_var_selecionadas, by = "variavel") %>% 
  mutate(ano = lubridate::ymd(ano, truncated = 2L))


IEPS_Completo_Municipio_SP = IEPS_Completo_2010_2022_selecionados_SP %>% 
  filter(nivel == "Município") %>% 
  pivot_longer(cols = cob_ab:pct_pop_masc,
               names_to = "variavel",
               values_to = "valor") %>% 
  inner_join(IEPS_var_selecionadas, by = "variavel") %>% 
  mutate(ano = lubridate::ymd(ano, truncated = 2L))


IEPS_Completo_Regiao_SP_vacinas = IEPS_Completo_Regiao_SP %>% 
  filter(str_detect(nome_dos_indicadores, "Cobertura Vacinal")) %>% 
  select(ano, regiao, no_regiao, nome_dos_indicadores, valor) %>% 
  mutate(valor = round(valor, 2),
         ano = as.character(ano) %>% str_replace(., "-01-01", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, "Cobertura Vacinal de ", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, " \\(\\%\\)", ""))


IEPS_Completo_Municipio_SP_vacinas = IEPS_Completo_Municipio_SP %>% 
  filter(str_detect(nome_dos_indicadores, "Cobertura Vacinal")) %>% 
  select(ano, regiao, nomemun, nome_dos_indicadores, valor) %>% 
  mutate(valor = round(valor, 2),
         ano = as.character(ano) %>% str_replace(., "-01-01", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, "Cobertura Vacinal de ", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, " \\(\\%\\)", ""))

Quais regiões tiveram cobertura vacinal mínima de 90% entre 2015 e 2022? Isto é, quais regiões foram menos afetadas pelo efeito 2016 e pela pandemia?

#Regiões -----
dados_vacinas = IEPS_Completo_Regiao_SP %>% 
  filter(str_detect(nome_dos_indicadores, "Cobertura Vacinal")) %>% 
  select(ano, regiao, no_regiao, nome_dos_indicadores, valor) %>% 
  mutate(valor = round(valor, 2),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, "Cobertura Vacinal de ", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, " \\(\\%\\)", "")) %>% 
  drop_na()

#Selecionar regiões com cobertura vacinal mínima de 90%
regioes_acima_90 = dados_vacinas %>% 
  filter(ano >= ymd("2015-01-01") & ano <= ymd("2022-01-01")) %>% 
  group_by(no_regiao, nome_dos_indicadores) %>% 
  summarise(min_coverage = min(valor)) %>%
  filter(min_coverage > 90)

#Municípios ------
dados_vacinas = IEPS_Completo_Municipio_SP %>% 
  filter(str_detect(nome_dos_indicadores, "Cobertura Vacinal")) %>% 
  select(ano, regiao, nomemun, nome_dos_indicadores, valor) %>% 
  mutate(valor = round(valor, 2),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, "Cobertura Vacinal de ", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, " \\(\\%\\)", "")) %>% 
  drop_na()

#Selecionar municipios com cobertura vacinal mínima de 90%
municipios_acima_90 = dados_vacinas %>% 
  filter(ano >= ymd("2015-01-01") & ano <= ymd("2022-01-01")) %>% 
  group_by(nomemun, nome_dos_indicadores) %>% 
  summarise(min_coverage = min(valor)) %>%
  filter(min_coverage > 90)

Quais foram as regiões que mais tiveram redução da cobertura vacinal por vacina?

# Calcular a variação anual e filtrar variações negativas
variacoes_negativas <- dados_vacinas %>%
    drop_na() %>% 
  arrange(no_regiao, nome_dos_indicadores, ano) %>%
  group_by(no_regiao, nome_dos_indicadores) %>%
  mutate(variacao_anual = valor - lag(valor)) %>%
  group_by(nome_dos_indicadores) %>% 
  arrange(nome_dos_indicadores, variacao_anual)  %>% 
  filter(variacao_anual<0)

#Top 5 regiões com maior queda da cobertura vacinal
variacoes_negativas_top5 = variacoes_negativas %>% 
  slice_min(n = 5, order_by = variacao_anual)
top5_reduz_2016 = variacoes_negativas %>% 
  filter(ano == "2016-01-01")  %>% 
  slice_min(n = 5, order_by = variacao_anual)
top5_reduz_2020 = variacoes_negativas %>% 
  filter(ano == "2020-01-01")  %>% 
  slice_min(n = 5, order_by = variacao_anual)

#Top 5 regiões com maior queda da cobertura de BCG
bcg_top5_reduz_2016 = variacoes_negativas %>% 
  filter(ano == "2016-01-01", 
         nome_dos_indicadores == "BCG")  %>% 
  slice_min(n = 5, order_by = variacao_anual)
bcg_top5_reduz_2020 = variacoes_negativas %>% 
  filter(ano == "2020-01-01", 
         nome_dos_indicadores == "BCG")  %>% 
  slice_min(n = 5, order_by = variacao_anual)

bcg_top5_reduz_2016
bcg_top5_reduz_2020

E os municípios?

# Calcular a variação anual e filtrar variações negativas
variacoes_negativas <- dados_vacinas %>%
    drop_na() %>% 
  arrange(nomemun, nome_dos_indicadores, ano) %>%
  group_by(nomemun, nome_dos_indicadores) %>%
  mutate(variacao_anual = valor - lag(valor)) %>%
  group_by(nome_dos_indicadores) %>% 
  arrange(nome_dos_indicadores, variacao_anual)  %>% 
  filter(variacao_anual<0)

#Top 5 regiões com maior queda da cobertura vacinal
variacoes_negativas_top5_mun = variacoes_negativas %>% 
  slice_min(n = 5, order_by = variacao_anual)
top5_reduz_2016_mun = variacoes_negativas %>% 
  filter(ano == "2016-01-01")  %>% 
  slice_min(n = 5, order_by = variacao_anual)
top5_reduz_2020_mun = variacoes_negativas %>% 
  filter(ano == "2020-01-01")  %>% 
  slice_min(n = 5, order_by = variacao_anual)

#Top 5 regiões com maior queda da cobertura de BCG
bcg_top5_reduz_2016_mun = variacoes_negativas %>% 
  filter(ano == "2016-01-01", 
         nome_dos_indicadores == "BCG")  %>% 
  slice_min(n = 5, order_by = variacao_anual)
bcg_top5_reduz_2020_mun = variacoes_negativas %>% 
  filter(ano == "2020-01-01", 
         nome_dos_indicadores == "BCG")  %>% 
  slice_min(n = 5, order_by = variacao_anual)

bcg_top5_reduz_2016_mun
bcg_top5_reduz_2020_mun

vacinas_estado_sp = IEPS_Completo_filtrado_SP %>% 
  filter(bloco == "Cobertura vacinal") %>% 
  bind_rows(brasil_vacinas) %>% 
  mutate(nivel = if_else(nivel == "Estado", "São Paulo", nivel),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, "Cobertura Vacinal de ", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, " \\(\\%\\)", ""),
         nivel = fct_rev(nivel)) %>% 
  ggplot() +
  aes(x = ano, y = valor, group = nivel) +
  scale_x_date(date_breaks = "1 year", date_labels = "%y") +
  scale_y_continuous(expand = expansion(add = 10)) +
  geom_rect(xmin = as.numeric(as.Date("2016-01-01")),
           xmax = as.numeric(as.Date("2020-01-01")),
            ymin = 0,
            ymax = 110,
            fill = "gray",
            alpha = 0.03) +
  geom_rect(xmin = as.numeric(as.Date("2020-01-01")),
               xmax = as.numeric(as.Date("2022-01-01")),
                ymin = 0,
                ymax = 110,
                fill = "gray40",
                alpha = 0.03) +
  geom_line(size = 2, mapping = aes(color = nivel, linetype = nivel)) +
  theme(text = element_text(size = 12),
        plot.title = element_text(size = 20),
        plot.subtitle = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.position = "top",
        strip.text.x = element_text(size = 12,
                                    margin = margin(5, 5, 15, 5, "pt"))) +
  labs(title = "Cobertura vacinal no estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "",
       y = "Cobertura vacinal (%)",
       x = "") +
  facet_wrap(~nome_dos_indicadores, scales = "free")

vacinas_estado_sp %>% 
  ggsave(file = here("Figuras", "Vacinas_Brasil_SP_2010_2022_linhas.png"), width = 12, height = 7)

vacinas_estado_sp

#Análise por microrregião

Densidade

IEPS_Completo_Regiao_SP_vacinas_ridge = IEPS_Completo_Regiao_SP_vacinas %>% 
  drop_na() %>% 
  ggplot() +
  aes(x = valor, y = fct_rev(ano), fill = ano) +
  geom_density_ridges() +
  scale_fill_manual(values = colorRampPalette(c("#65ADC2", "purple"))(15)) +
  theme(legend.position = "none",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        strip.text.x = element_text(size = 15, face = "bold"),
        panel.spacing=unit(2, "lines")) +
  labs(title = "Cobertura vacinal",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       y = "",
       x = "Cobertura vacinal (%)") +
    facet_wrap(~nome_dos_indicadores, scales = "free_x", nrow = 2)
  
IEPS_Completo_Regiao_SP_vacinas_ridge %>% 
    ggsave(file = here("Figuras", "vacinas_regioes_estado_sp_2010_2022_densidade.png"), width = 15, height = 10)

IEPS_Completo_Regiao_SP_vacinas_ridge

Boxplots

 IEPS_Completo_Regiao_SP_vacinas_boxplot = IEPS_Completo_Regiao_SP_vacinas %>% 
  drop_na() %>% 
  ggplot() +
  aes(x = valor, y = fct_rev(ano), fill = ano) +
  geom_boxplot(outliers = F) +
  geom_jitter(aes(label = no_regiao, color = valor),
              alpha = 0.2,
              na.rm = T) +
  scale_fill_manual(values = colorRampPalette(c("#65ADC2", "purple"))(15)) +
  facet_wrap(vars(nome_dos_indicadores))  +
    theme(legend.position = "none",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        strip.text.x = element_text(size = 15, face = "bold"),
        panel.spacing=unit(2, "lines")) +
  labs(title = "Cobertura vacinal",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       y = "",
       x = "Cobertura vacinal (%)") +
    facet_wrap(~nome_dos_indicadores, scales = "free_x", nrow = 2)
IEPS_Completo_Regiao_SP_vacinas_boxplot %>% 
    ggsave(file = here("Figuras", "vacinas_regioes_estado_sp_2010_2022_boxplot.png"), width = 15, height = 12)

IEPS_Completo_Regiao_SP_vacinas_boxplot  

IEPS_Completo_Regiao_SP_vacinas_BCG_boxplot = IEPS_Completo_Regiao_SP_vacinas %>% 
  filter(nome_dos_indicadores == "BCG") %>% 
  drop_na() %>% 
  ggplot() +
  aes(x = fct_rev(ano), y = valor, fill = ano) +
  geom_boxplot(outliers = F) +
  geom_jitter(aes(label = no_regiao, color = valor),
              alpha = 0.2,
              na.rm = T) +
  scale_fill_manual(values = colorRampPalette(c("#65ADC2", "purple"))(15)) +
    theme(legend.position = "none",
        text = element_text(size = 12, color = "black"),
        plot.title = element_text(size = 20),
        plot.subtitle = element_text(size = 15),
        plot.caption = element_text(size = 10),
        plot.margin = unit(c(1,1,1,1), "cm"),
        strip.text.x = element_text(size = 10, face = "bold"),
        panel.spacing=unit(2, "lines")) +
  labs(title = "Cobertura vacinal - BCG",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       y = "",
       x = "Cobertura vacinal (%)")  

ggplotly(IEPS_Completo_Regiao_SP_vacinas_BCG_boxplot)

Gráfico de linhas

regioes_acima_90_vacina_linhas = IEPS_Completo_Regiao_SP %>% 
  filter(str_detect(nome_dos_indicadores, "Cobertura Vacinal"),
         no_regiao %in% unique(regioes_acima_90$no_regiao)) %>% 
  select(ano, regiao, no_regiao, nome_dos_indicadores, valor) %>% 
  mutate(valor = round(valor, 2),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, "Cobertura Vacinal de ", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, " \\(\\%\\)", "")) %>% 
  drop_na() %>% 
  ggplot() +
  aes(x = ano, y = valor, group = no_regiao) +
  geom_line(aes(color = no_regiao),
            size = 2) +
  scale_x_date(date_breaks = "1 year", date_labels = "%y") +
    theme(legend.position = "top",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal no estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Microrregião",
       y = "Cobertura vacinal (%)",
       x = "") +
  facet_wrap(~nome_dos_indicadores, scales = "free") +
  scale_shape_manual(values = c(16)) +
  scale_color_cosmic()

regioes_acima_90_vacina_linhas
regioes_acima_90_vacina_linhas %>% 
    ggsave(file = here("Figuras", "vacinas_acima_90_regioes_estado_sp_2010_2022_linhas.png"), width = 15, height = 10)

bcg_queda_2016_linhas = IEPS_Completo_Regiao_SP %>% 
  filter(str_detect(nome_dos_indicadores, "Cobertura Vacinal"),
         str_detect(nome_dos_indicadores, "BCG"),
         no_regiao %in% bcg_top5_reduz_2016$no_regiao) %>% 
  select(ano, regiao, no_regiao, nome_dos_indicadores, valor) %>% 
  mutate(valor = round(valor, 2),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, "Cobertura Vacinal de ", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, " \\(\\%\\)", "")) %>% 
  drop_na() %>% 
  ggplot() +
  aes(x = ano, y = valor, group = no_regiao) +
  geom_line(aes(color = no_regiao),
            size = 2) +
  scale_x_date(date_breaks = "1 year", date_labels = "%y") +
    theme(legend.position = "right",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal da BCG no estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Microrregião",
       y = "Cobertura vacinal (%)",
       x = "") +
  scale_color_cosmic()
bcg_queda_2016_linhas

bcg_queda_2016_linhas %>% 
    ggsave(file = here("Figuras", "vacinas_queda_2016_regioes_estado_sp_2010_2022_linhas.png"), width = 10, height = 5)

2020

bcg_queda_2020_linhas = IEPS_Completo_Regiao_SP %>% 
  filter(str_detect(nome_dos_indicadores, "Cobertura Vacinal"),
         str_detect(nome_dos_indicadores, "BCG"),
         no_regiao %in% bcg_top5_reduz_2020$no_regiao) %>% 
  select(ano, regiao, no_regiao, nome_dos_indicadores, valor) %>% 
  mutate(valor = round(valor, 2),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, "Cobertura Vacinal de ", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, " \\(\\%\\)", "")) %>% 
  drop_na() %>% 
  ggplot() +
  aes(x = ano, y = valor, group = no_regiao) +
  geom_line(aes(color = no_regiao),
            size = 2) +
  scale_x_date(date_breaks = "1 year", date_labels = "%y") +
    theme(legend.position = "right",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal da BCG no estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Microrregião",
       y = "Cobertura vacinal (%)",
       x = "") +
  scale_color_cosmic()
bcg_queda_2020_linhas

bcg_queda_2020_linhas %>% 
    ggsave(file = here("Figuras", "vacinas_queda_2020_regioes_estado_sp_2010_2022_linhas.png"), width = 10, height = 5)

#Análise por municipio

Densidade

Boxplots

 IEPS_Completo_Municipio_SP_vacinas_boxplot = IEPS_Completo_Municipio_SP %>% 
  filter(str_detect(variavel, "cob_vac_")) %>% 
  drop_na() %>% 
    mutate(ano = as.character(ano)) %>% 
  ggplot() +
  aes(x = valor, y = fct_rev(ano), fill = ano) +
  geom_jitter(aes(label = no_regiao, color = valor),
              alpha = 0.2,
              na.rm = T) +
  geom_boxplot(outliers = F) +
  scale_fill_manual(values = colorRampPalette(c("#65ADC2", "purple"))(15)) +
  theme(legend.position = "none",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        strip.text.x = element_text(size = 15, face = "bold"),
        panel.spacing=unit(2, "lines")) +
  labs(title = "Cobertura vacinal, por municipio",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       y = "",
       x = "Cobertura vacinal (%)") +
    facet_wrap(~nome_dos_indicadores, scales = "free_x", nrow = 2)


IEPS_Completo_Municipio_SP_vacinas_boxplot %>% 
    ggsave(file = here("Figuras", "Vacinas_Municipios_SP_2010_2022_boxplot.png"), width = 20, height = 12)

IEPS_Completo_Municipio_SP_vacinas_boxplot  

Boxplot interativo

IEPS_Completo_Regiao_SP_vacinas_BCG_boxplot = IEPS_Completo_Regiao_SP_vacinas %>% 
  filter(nome_dos_indicadores == "BCG") %>% 
  drop_na() %>% 
  ggplot() +
  aes(x = fct_rev(ano), y = valor, fill = ano) +
  geom_boxplot(outliers = F) +
  geom_jitter(aes(label = no_regiao, color = valor),
              alpha = 0.2,
              na.rm = T) +
  scale_fill_manual(values = colorRampPalette(c("#65ADC2", "purple"))(15)) +
    theme(legend.position = "none",
        text = element_text(size = 12, color = "black"),
        plot.title = element_text(size = 20),
        plot.subtitle = element_text(size = 15),
        plot.caption = element_text(size = 10),
        plot.margin = unit(c(1,1,1,1), "cm"),
        strip.text.x = element_text(size = 10, face = "bold"),
        panel.spacing=unit(2, "lines")) +
  labs(title = "Cobertura vacinal - BCG",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       y = "",
       x = "Cobertura vacinal (%)")  

ggplotly(IEPS_Completo_Regiao_SP_vacinas_BCG_boxplot)

Gráfico de linhas

regioes_acima_90_vacina_linhas = IEPS_Completo_Regiao_SP %>% 
  filter(str_detect(nome_dos_indicadores, "Cobertura Vacinal"),
         no_regiao %in% unique(regioes_acima_90$no_regiao)) %>% 
  select(ano, regiao, no_regiao, nome_dos_indicadores, valor) %>% 
  mutate(valor = round(valor, 2),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, "Cobertura Vacinal de ", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, " \\(\\%\\)", "")) %>% 
  drop_na() %>% 
  ggplot() +
  aes(x = ano, y = valor, group = no_regiao) +
  geom_line(aes(color = no_regiao),
            size = 2) +
  scale_x_date(date_breaks = "1 year", date_labels = "%y") +
    theme(legend.position = "top",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal no estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Microrregião",
       y = "Cobertura vacinal (%)",
       x = "") +
  facet_wrap(~nome_dos_indicadores, scales = "free") +
  scale_shape_manual(values = c(16)) +
  scale_color_cosmic()

regioes_acima_90_vacina_linhas
regioes_acima_90_vacina_linhas %>% 
    ggsave(file = here("Figuras", "vacinas_acima_90_regioes_estado_sp_2010_2022_linhas.png"), width = 15, height = 10)

bcg_queda_2016_linhas = IEPS_Completo_Regiao_SP %>% 
  filter(str_detect(nome_dos_indicadores, "Cobertura Vacinal"),
         str_detect(nome_dos_indicadores, "BCG"),
         no_regiao %in% bcg_top5_reduz_2016$no_regiao) %>% 
  select(ano, regiao, no_regiao, nome_dos_indicadores, valor) %>% 
  mutate(valor = round(valor, 2),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, "Cobertura Vacinal de ", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, " \\(\\%\\)", "")) %>% 
  drop_na() %>% 
  ggplot() +
  aes(x = ano, y = valor, group = no_regiao) +
  geom_line(aes(color = no_regiao),
            size = 2) +
  scale_x_date(date_breaks = "1 year", date_labels = "%y") +
    theme(legend.position = "right",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal da BCG no estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Microrregião",
       y = "Cobertura vacinal (%)",
       x = "") +
  scale_color_cosmic()
bcg_queda_2016_linhas

bcg_queda_2016_linhas %>% 
    ggsave(file = here("Figuras", "vacinas_queda_2016_regioes_estado_sp_2010_2022_linhas.png"), width = 10, height = 5)

2020

bcg_queda_2020_linhas = IEPS_Completo_Regiao_SP %>% 
  filter(str_detect(nome_dos_indicadores, "Cobertura Vacinal"),
         str_detect(nome_dos_indicadores, "BCG"),
         no_regiao %in% bcg_top5_reduz_2020$no_regiao) %>% 
  select(ano, regiao, no_regiao, nome_dos_indicadores, valor) %>% 
  mutate(valor = round(valor, 2),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, "Cobertura Vacinal de ", ""),
         nome_dos_indicadores = str_replace(nome_dos_indicadores, " \\(\\%\\)", "")) %>% 
  drop_na() %>% 
  ggplot() +
  aes(x = ano, y = valor, group = no_regiao) +
  geom_line(aes(color = no_regiao),
            size = 2) +
  scale_x_date(date_breaks = "1 year", date_labels = "%y") +
    theme(legend.position = "right",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal da BCG no estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Microrregião",
       y = "Cobertura vacinal (%)",
       x = "") +
  scale_color_cosmic()
bcg_queda_2020_linhas

bcg_queda_2020_linhas %>% 
    ggsave(file = here("Figuras", "vacinas_queda_2020_regioes_estado_sp_2010_2022_linhas.png"), width = 10, height = 5)

PCA

Microrregião

library(tidyverse)
library(ggrepel) #load before factoextra
library(factoextra)
library(caret)
library(stats)
library(ggfortify)
library(ggplot2)
library(corrplot)
library(cowplot)
library(ggcorrplot)
library(ComplexHeatmap)
library(circlize)
library(ggExtra)
library(FactoMineR)

Microrregião: Selecionados

#Imputs
ano_selec = "2020"
ano_selec_ymd = ymd("2020-01-01")
row_names = "no_regiao"
nivel_selec = "Região"

#
pca_dados = IEPS_Completo_2010_2022_selecionados_SP %>% 
  filter(nivel == nivel_selec,
         ano == ano_selec) %>%
  column_to_rownames(row_names) %>% 
  select(cob_ab:cob_vac_penta, 
         tx_mort_aj_cens:pct_desp_recp_saude_mun, 
         desp_tot_saude_pc_mun_def,
         desp_recp_saude_pc_mun_def, 
         ideb_5ano:pop, 
         pct_pop_fem:pct_pop_masc)

pca_dados_matrix = pca_dados %>% 
  as.matrix()

dim(pca_dados_matrix)

pca_dados_ann = IEPS_Completo_2010_2022_selecionados_SP %>% 
  filter(nivel == nivel_selec,
         ano == ano_selec) %>% 
  select(row_names, pop) %>% 
  mutate(porte_municipio = case_when(pop <= 20000 ~ "Pequeno Porte I",
                                     pop <= 50000 ~ "Pequeno Porte II",
                                     pop <= 100000 ~ "Médio Porte",
                                     pop <= 900000 ~ "Grande Porte",
                                     TRUE ~ "Metrópole"),
         ID = row_names) %>% 
  select(-pop) %>% 
    column_to_rownames(row_names)

# Delete columns with near 0 variance
nearZeroVarCols <- nearZeroVar(pca_dados_matrix, saveMetrics = TRUE)
pca_dados_matrix_pca <- pca_dados_matrix[, !nearZeroVarCols$nzv]
pca_res <- prcomp(pca_dados_matrix_pca, scale. = T)

#Selecionar variáveis de interesse
corr_matrix_selec = corr_matrix %>% 
  as.data.frame() %>% 
  rownames_to_column("id") %>% 
  filter(str_detect(id, "cob_vac")) %>% 
  column_to_rownames("id") %>% 
  as.matrix()

# Correlation plot ----
corr_matrix = cor(pca_dados_matrix_pca %>% scale()) 
var <- get_pca_var(pca_res)
p.mat = cor_pmat(pca_dados_matrix_pca %>% scale())

p.mat_2 = cor_pmat(corr_matrix_selec %>% scale()) %>% 
  as.data.frame() %>% 
  rownames_to_column("id") %>% 
  inner_join(corr_matrix_selec %>% as.data.frame() %>% rownames_to_column("id") %>% select("id"), by = "id") %>% 
  column_to_rownames("id") %>% 
  as.matrix()



#Selecionados
png(file = here("Figuras", paste0("PCA_", nivel_selec ,"_", ano_selec, "_Selecionados_Corrplot.png")), 
    width = 15, height = 5,
    units = "in",
    res = 300)

corrplot(round(corr_matrix_selec, 2),
         is.corr = TRUE,
         tl.col = "black",
         title = "Cobertura vacinal e indicadores sociais, economicos, demográficos e de saúde",
         cex.main = 2,
         method = "color",
         p.mat = p.mat_2,
         tl.srt = 45,
         sig.level = c(0.05),
         mar=c(1,0,4,0),
         addCoef.col = 'black',
         number.cex = 0.6,
         insig = "blank",
         col = colorRampPalette(c("#65ADC2", "white", "#E84646"))(11))

dev.off()

#Todos
corrplot = ggcorrplot(corr_matrix,
                      hc.order = TRUE,
                      method = "circle",
                      type = "full",
                      ggtheme = NULL,
                      outline.col = "white",
                      p.mat = p.mat,
         insig = "blank",

                      sig.level = 0.05,
                      colors = c("#65ADC2", "white", "#E84646")) +
  theme(axis.text.x = element_text(angle = 45, size = 10),
        axis.text.y = element_text(size = 10)) +
  labs(title = "Cobertura vacinal e indicadores sociais, economicos, demográficos e de saúde") +
  theme(plot.title = element_text(size = 15, face = "bold"))

corrplot %>% 
  ggsave(file = here("Figuras", paste0("PCA_", nivel_selec ,"_", ano_selec, "_TodosIndicadores_Corrplot.png")), width = 15, height = 10)

#PCA ----
data.pca = prcomp(corr_matrix) #PCA
summary(corr_matrix) #Retornar PCs

#Scree plot ----
scree_plot = fviz_eig(data.pca, 
                      addlabels = TRUE,
                      ylim = c(0, 70)) +
  geom_col(color = "#00AFBB", fill = "#00AFBB") 

scree_plot

scree_plot %>% 
  ggsave(file = here("Figuras", paste0("PCA_", nivel_selec ,"_", ano_selec, "_Screeplot.png")), 
         width = 5,
         height = 3)

# Loadings plot ----
options(ggrepel.max.overlaps = Inf)

circle_contrib= fviz_pca_var(pca_res, col.var = "cos2",
                              gradient.cols = c("#65ADC2", "black"),
                              select.var= list(cos2 = 30), 
                              repel = T, 
                              labelsize = 4,
                              col.circle = NA) 

circle_contrib

circle_contrib %>% 
  ggsave(file = here("Figuras", paste0("PCA_", nivel_selec ,"_", ano_selec, "_Loadings.png")), 
         width = 5, 
         height = 5)

# Cluster ------
#Determine the number of clusters
pca_scores <- data.frame(pca_res$x[, 1:2])
fviz_nbclust(pca_scores,  
                     FUNcluster = kmeans,
                     method = "wss")


set.seed(666) # Set seed for randomization
cluster_model <- kmeans(pca_res$x[, 1:2], centers = 2)  #Set the number of clusters

pca_dados$cluster <- as.factor(cluster_model$cluster)
pca_dados = pca_dados %>% 
  rownames_to_column("id")

# Plot ------
pca_plot_knn_cluster = autoplot(pca_res, 
        data = pca_dados, 
        colour = 'cluster')  +
  stat_ellipse(aes(color = cluster, 
                   fill = cluster), 
               geom = "polygon", 
               alpha = 0.1, 
               linetype = 1,
               size = 0.3,
               type = "t") +
  geom_point(aes(fill = cluster,
                 size = pop),
             colour="black",
             pch=21) +
  labs(title = paste0(nivel_selec, "paulista, ", ano_selec)) + 
  theme(text = element_text(color = "black")) +
  geom_label_repel(aes(label = id,
                      segment.colour= "black"),
                  box.padding = 0.3,
                  size = 3)  +
  scale_size_continuous(range = c(2,8))
pca_plot_knn_cluster_marginal = ggMarginal(
  pca_plot_knn_cluster,
  groupFill = T,
  groupColour = T)
pca_plot_knn_cluster_marginal

pca_plot_knn_cluster_marginal %>% ggsave(file = here("Figuras", paste0("PCA_", nivel_selec ,"_", ano_selec, "_Clusters.png")), width = 10, height = 6)

ggthemr("fresh", spacing = 1)
swatch()

#Polio ----
pca_plot_knn_cluster_gradient = autoplot(pca_res, 
        data = pca_dados)  +
  geom_point(aes(fill = cob_vac_polio,
                 size = pop),
             colour="black",
             pch=21) +
  labs(title = paste0(nivel_selec, "paulista, ", ano_selec)) + 
  theme(text = element_text(color = "black"),
        legend.position = "right")  +
  scale_fill_gradient(low = "#E84646", high = "#65ADC2") +
    scale_size_continuous(range = c(2,8))


pca_plot_knn_cluster_gradient

pca_plot_knn_cluster_gradient %>% ggsave(file =  here("Figuras", paste0("PCA_", nivel_selec ,"_", ano_selec,"_Polio.png")), width = 10, height = 6)

#BCG -----

pca_plot_knn_cluster_gradient = autoplot(pca_res, 
        data = pca_dados)  +
  geom_point(aes(fill = cob_vac_bcg,
                 size = pop),
             colour="black",
             pch=21) +
  labs(title = paste0(nivel_selec, "paulista, ", ano_selec)) + 
  theme(text = element_text(color = "black"),
        legend.position = "right")  +
  scale_fill_gradient(low = "#E84646", high = "#65ADC2") +
    scale_size_continuous(range = c(2,8))


pca_plot_knn_cluster_gradient

pca_plot_knn_cluster_gradient %>% ggsave(file =  here("Figuras", paste0("PCA_", nivel_selec ,"_", ano_selec,"_BCG.png")), width = 10, height = 6)

"#111111" "#65ADC2" "#233B43" "#E84646" "#C29365" "#362C21" "#316675" "#168E7F" "#109B37"
res.km <- kmeans(var$coord, centers = 3, nstart = 25)
grp <- as.factor(res.km$cluster)

pca_plot_knn_cluster_biplot = fviz_pca_biplot(pca_res, 
                col.var = "black",
                geom.var = c("point", "text"),
                fill.ind = pca_dados$cluster,
                col.ind = "black",
                pointshape = 21,
                palette = c("#65ADC2", "#233B43", "#E84646"),
                label = "var",
                select.var= list(cos2 = 20), 
                pointsize = 3,
                addEllipses = T,
                repel = TRUE,
                legend.title = "Cluster") 

pca_plot_knn_cluster_biplot

pca_plot_knn_cluster_biplot %>% ggsave(file = here("Figuras", paste0("PCA_", nivel_selec ,"_", ano_selec,"_Biplot.png")), width = 10, height = 6)

Microrregião: Dispersão

pop_regiao = IEPS_Completo_2010_2022_selecionados_SP %>%
  filter(nivel == nivel_selec) %>% 
    mutate(porte_regiao = case_when(pop <= 20000 ~ "Pequeno Porte I",
                                     pop <= 50000 ~ "Pequeno Porte II",
                                     pop <= 100000 ~ "Médio Porte",
                                     pop <= 900000 ~ "Grande Porte",
                                     TRUE ~ "Metrópole")) %>% 
  select(row_names, porte_municipio)


dados_vacinas = IEPS_Completo_2010_2022_selecionados_SP %>% 
  filter(nivel == nivel_selec) %>%
  select(row_names, ano, cob_ab:cob_vac_penta, 
         tx_mort_aj_cens:pct_desp_recp_saude_mun, 
         desp_tot_saude_pc_mun_def,
         desp_recp_saude_pc_mun_def, 
         ideb_5ano:pop, 
         pct_pop_fem:pct_pop_masc) %>% 
  inner_join(pop_regiao, by = row_names) %>% 
  distinct()

Correlação ideb e vacina

# library(esquisse)
# library(ggpmisc)
cor_vac_ideb = dados_vacinas %>% 
  filter(!ano == ymd("2022-01-01")) %>% 
  select(c(row_names, "ideb_5ano", "ideb_9ano", "cob_vac_bcg", "ano", "porte_municipio", "pop")) %>% 
  ggplot(.) +
  aes(x = ideb_5ano, y = cob_vac_bcg) +
  # geom_point(color = "black", 
  #            shape = "circle") +
  geom_point(aes(fill = porte_municipio, 
                 colour = porte_municipio,
                 size = pop), 
             shape = "circle") +
  geom_smooth(method = "lm",
              alpha = 0.1,
              se = T) +
  scale_y_continuous(expand = expansion(add = 20)) +
      scale_size_continuous(range = c(2,8)) +
  stat_correlation(method = "pearson",
                   size = 3.5,
                   label.x = 5, 
                   label.y = 120,
                   mapping = use_label(c("R", "R2", "P", "n"))) +
  theme(legend.position = "right",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal da BCG e IDEB 5o ano, microrregiões do estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Microrregião",
       y = "Cobertura vacinal (%)",
       x = "IDEB") +
  facet_wrap(~ano, scales = "free", nrow = 2)

cor_vac_ideb

cor_vac_ideb %>% 
  ggsave(file = here("Figuras", here("Figuras", paste0("PCA_", nivel_selec ,"_", ano_selec,"_2010_2020_BCG_vs_IDEB.png")) 
         width = 20,
         height = 10)

IDHM e vacina

cor_vac_idhm = dados_vacinas %>% 
  filter(!ano == ymd("2022-01-01")) %>% 
  select(c("no_regiao", "idhm", "cob_vac_bcg", "ano", "porte_municipio", "pop")) %>% 
  ggplot(.) +
  aes(x = idhm, y = cob_vac_bcg) +
  # geom_point(color = "black", 
  #            shape = "circle") +
  geom_point(aes(fill = porte_municipio, 
                 colour = porte_municipio,
                 size = pop), 
             shape = "circle") +
  geom_smooth(method = "lm",
              alpha = 0.1,
              se = T) +
  scale_y_continuous(expand = expansion(add = 20)) +
      scale_size_continuous(range = c(2,8)) +
  stat_correlation(method = "pearson",
                   size = 3.5,
                   label.x = 5, 
                   label.y = 120,
                   mapping = use_label(c("R", "R2", "P", "n"))) +
  theme(legend.position = "right",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal da BCG e IDH municipal, microrregiões do estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Microrregião",
       y = "Cobertura vacinal (%)",
       x = "IDHM") +
  facet_wrap(~ano, scales = "free", nrow = 2)

cor_vac_idhm

cor_vac_idhm %>% 
  ggsave(file = here("Figuras", "Microrregioes_2010_2020_BCG_vs_IDHM.png"), 
         width = 20, height = 10)

Taxa de mortes evitaveis e vacina

variavel_interesse = "tx_mort_evit_aj_cens"

cor_vac_txmortes = dados_vacinas %>% 
  filter(!ano == ymd("2022-01-01")) %>% 
  select(c("no_regiao", variavel_interesse, "cob_vac_bcg", "ano", "porte_municipio", "pop")) %>% 
  ggplot(.) +
  aes(x = variavel_interesse, y = cob_vac_bcg) +
  # geom_point(color = "black", 
  #            shape = "circle") +
  geom_point(aes(fill = porte_municipio, 
                 colour = porte_municipio,
                 size = pop), 
             shape = "circle") +
  geom_smooth(method = "lm",
              alpha = 0.1,
              se = T) +
  scale_y_continuous(expand = expansion(add = 20)) +
      scale_size_continuous(range = c(2,8)) +
  stat_correlation(method = "pearson",
                   size = 3.5,
                   label.x = 5, 
                   label.y = 120,
                   mapping = use_label(c("R", "R2", "P", "n"))) +
  theme(legend.position = "right",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal da BCG e Taxa de mortes evitaveis (CENSO), microrregiões do estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Microrregião",
       y = "Taxa de mortes evitaveis (CENSO)",
       x = "IDHM") +
  facet_wrap(~ano, scales = "free", nrow = 2)

cor_vac_txmortes

cor_vac_txmortes %>% 
  ggsave(file = here("Figuras", "Microrregioes_2010_2020_BCG_vs_Taxa_Mortes_Evitaveis.png"), 
         width = 20, height = 10)

Município

library(tidyverse)
library(ggrepel) #load before factoextra
library(factoextra)
library(caret)
library(stats)
library(ggfortify)
library(ggplot2)
library(corrplot)
library(cowplot)
library(ggcorrplot)
library(ComplexHeatmap)
library(circlize)
library(ggExtra)
library(FactoMineR)

Muncípios: Selecionados

ano_selec = "2020-01-01"
nivel_selec = "Município"
row_names = "nomemun"

pca_municipio = IEPS_Completo_2010_2022_selecionados_SP %>% 
  filter(nivel == "Município",
         ano == ymd(ano_selec)) %>%
  column_to_rownames("nomemun") %>% 
  select(cob_ab:cob_vac_penta, 
         tx_mort_aj_cens:pct_desp_recp_saude_mun, 
         desp_tot_saude_pc_mun_def,
         desp_recp_saude_pc_mun_def, 
         ideb_5ano:pop, 
         pct_pop_fem:pct_pop_masc)

pca_municipio_matrix = pca_municipio %>% 
  as.matrix()

dim(pca_municipio_matrix)

pca_municipio_ann = IEPS_Completo_2010_2022_selecionados_SP %>% 
  filter(nivel == "Município",
         ano == ymd(ano_selec)) %>% 
  select(nomemun, pop) %>% 
  mutate(porte_municipio = case_when(pop <= 20000 ~ "Pequeno Porte I",
                                     pop <= 50000 ~ "Pequeno Porte II",
                                     pop <= 100000 ~ "Médio Porte",
                                     pop <= 900000 ~ "Grande Porte",
                                     TRUE ~ "Metrópole"),
         ID = nomemun) %>% 
  select(-pop) %>% 
    column_to_rownames("nomemun")

pca_municipio %>% 
  ggplot()

# Delete columns with near 0 variance
nearZeroVarCols <- nearZeroVar(pca_municipio_matrix, saveMetrics = TRUE)
pca_municipio_matrix_pca <- pca_municipio_matrix[, !nearZeroVarCols$nzv]
pca_res <- prcomp(pca_municipio_matrix_pca, scale. = T)

# Correlation plot ----
corr_matrix = cor(pca_municipio_matrix_pca %>% scale()) 
var <- get_pca_var(pca_res)

corr_matrix_selec = corr_matrix %>% 
  as.data.frame() %>% 
  rownames_to_column("id") %>% 
  filter(str_detect(id, "cob_vac")) %>% 
  column_to_rownames("id") %>% 
  as.matrix()


p.mat = cor_pmat(pca_municipio_matrix_pca %>% scale())

p.mat_2 = cor_pmat(corr_matrix_selec %>% scale()) %>% 
  as.data.frame() %>% 
  rownames_to_column("id") %>% 
  inner_join(corr_matrix_selec %>% as.data.frame() %>% rownames_to_column("id") %>% select("id"), by = "id") %>% 
  column_to_rownames("id") %>% 
  as.matrix()



#Selecionados
png(file = here("Figuras", paste0("PCA_", nivel_selec ,"_", ano_selec, "_Selecionados_Corrplot.png")), 
    width = 15, height = 5,
    units = "in",
    res = 300)

corrplot(round(corr_matrix_selec, 2),
         is.corr = TRUE,
         tl.col = "black",
         title = "Cobertura vacinal e indicadores sociais, economicos, demográficos e de saúde",
         cex.main = 2,
         method = "color",
         p.mat = p.mat_2,
         tl.srt = 45,
         sig.level = c(0.05),
         mar=c(1,0,4,0),
         addCoef.col = 'black',
         number.cex = 0.6,
         insig = "blank",
         col = colorRampPalette(c("#65ADC2", "white", "#E84646"))(11))

dev.off()

#Todos os indicadores

corrplot = ggcorrplot(corr_matrix, hc.order = TRUE,
           method = "circle",
           type = "full",
           ggtheme = NULL,
           outline.col = "white", 
           p.mat = p.mat,
           sig.level = 0.05, 
           colors = c("#65ADC2", "white", "#E84646")) +
  theme(axis.text.x = element_text(angle = 45, size = 10),
        axis.text.y = element_text(size = 10)) +
  labs(title = "Cobertura vacinal e indicadores sociais, economicos, demográficos e de saúde") +
  theme(plot.title = element_text(size = 15, face = "bold"))

corrplot

corrplot %>% 
  ggsave(file = here("Figuras", "PCA_Municípios_2020_Corrplot.png"), width = 12, height = 10)

#PCA ----
data.pca = prcomp(corr_matrix) #PCA
summary(corr_matrix) #Retornar PCs

#Scree plot ----
scree_plot = fviz_eig(data.pca, 
                      addlabels = TRUE,
                      ylim = c(0, 70)) +
  geom_col(color = "#00AFBB", fill = "#00AFBB") 

scree_plot

scree_plot %>% 
  ggsave(file = here("Figuras", "PCA_Municípios_2020_Selecionados_Screeplot.png"), 
         width = 5,
         height = 3)

# Loadings plot ----
options(ggrepel.max.overlaps = Inf)

circle_contrib= fviz_pca_var(pca_res, col.var = "cos2",
                              gradient.cols = c("#65ADC2", "black"),
                              select.var= list(cos2 = 30), 
                              repel = T, 
                              labelsize = 4,
                              col.circle = NA) 

circle_contrib

circle_contrib %>% 
  ggsave(file = here("Figuras", "PCA_Municípios_2020_Selecionados_Loadings.png"), 
         width = 5, 
         height = 5)

# Cluster ------
#Determine the number of clusters
pca_scores <- data.frame(pca_res$x[, 1:2])
fviz_nbclust(pca_scores,  
                     FUNcluster = kmeans,
                     method = "wss")


set.seed(666) # Set seed for randomization
cluster_model <- kmeans(pca_res$x[, 1:2], centers = 3)  #Set the number of clusters

pca_municipio$cluster <- as.factor(cluster_model$cluster)
pca_municipio = pca_municipio %>% 
  rownames_to_column("id")

# Plot ------
pca_plot_knn_cluster = autoplot(pca_res, 
        data = pca_municipio, 
        colour = 'cluster')  +
  stat_ellipse(aes(color = cluster, 
                   fill = cluster), 
               geom = "polygon", 
               alpha = 0.1, 
               linetype = 1,
               size = 0.3,
               type = "t") +
  geom_point(aes(fill = cluster,
                 size = pop),
             colour="black",
             pch=21) +
  labs(title = "Municípios paulistas, 2020") + 
  theme(text = element_text(color = "black")) +
  # geom_label_repel(aes(label = id,
  #                     segment.colour= "black"),
  #                 box.padding = 0.3,
  #                 size = 3)  +
  scale_size_continuous(range = c(1,8))
pca_plot_knn_cluster_marginal = ggMarginal(
  pca_plot_knn_cluster,
  groupFill = T,
  groupColour = T)
pca_plot_knn_cluster_marginal

pca_plot_knn_cluster_marginal %>% ggsave(file = here("Figuras", "PCA_Municípios_2020_Selecionados_Clusters.png"), width = 10, height = 6)

ggthemr("fresh", spacing = 1)
swatch()

pca_plot_knn_cluster_gradient = autoplot(pca_res, 
        data = pca_municipio)  +
  geom_point(aes(fill = cob_vac_polio,
                 size = pop),
             colour="black",
             pch=21) +
  labs(title = "Microrregiões paulistas, 2020") + 
  theme(text = element_text(color = "black"),
        legend.position = "right")  +
  scale_fill_gradient(low = "#E84646", high = "#65ADC2") +
    scale_size_continuous(range = c(4,8))


pca_plot_knn_cluster_gradient

pca_plot_knn_cluster_gradient %>% ggsave(file = here("Figuras", "PCA_Municipios_2020_Selecionados_Polio.png"), width = 10, height = 6)

ggthemr("fresh", spacing = 1)
swatch()

pca_municipio_2 = pca_municipio %>% 
  bind_cols(pca_municipio_ann)
pca_plot_knn_cluster_gradient = autoplot(pca_res, 
        data = pca_municipio_2)  +
  geom_point(aes(fill = porte_municipio,
                 size = pop),
             colour="black",
             pch=21) +
  labs(title = "Microrregiões paulistas, 2020") + 
  theme(text = element_text(color = "black"),
        legend.position = "right")  +
    scale_size_continuous(range = c(4,8))


pca_plot_knn_cluster_gradient

pca_plot_knn_cluster_gradient %>% ggsave(file = here("Figuras", "PCA_Municipios_2020_Selecionados_PorteMunicpio.png"), width = 10, height = 6)

"#111111" "#65ADC2" "#233B43" "#E84646" "#C29365" "#362C21" "#316675" "#168E7F" "#109B37"
res.km <- kmeans(var$coord, centers = 3, nstart = 25)
grp <- as.factor(res.km$cluster)

pca_plot_knn_cluster_biplot = fviz_pca_biplot(pca_res, 
                col.var = "black",
                geom.var = c("point", "text"),
                fill.ind = pca_municipio$cluster,
                col.ind = "black",
                pointshape = 21,
                palette = c("#65ADC2", "#233B43", "#E84646"),
                label = "var",
                select.var= list(cos2 = 20), 
                pointsize = 3,
                addEllipses = T,
                repel = TRUE,
                legend.title = "Cluster") 

pca_plot_knn_cluster_biplot

pca_plot_knn_cluster_biplot %>% ggsave(file = here("Figuras", "PCA_Municipios_2020_Selecionados_Biplot.png"), width = 10, height = 6)

Correlação: Dispersão

pop_municipio = IEPS_Completo_2010_2022_selecionados_SP %>%
  filter(nivel == "Município") %>% 
    mutate(porte_municipio = case_when(pop <= 20000 ~ "Pequeno Porte I",
                                     pop <= 50000 ~ "Pequeno Porte II",
                                     pop <= 100000 ~ "Médio Porte",
                                     pop <= 900000 ~ "Grande Porte",
                                     TRUE ~ "Metrópole")) %>% 
  select(nomemun, porte_municipio)


dados_vacinas = IEPS_Completo_2010_2022_selecionados_SP %>% 
  filter(nivel == "Município") %>%
  select(nomemun, ano, cob_ab:cob_vac_penta, 
         tx_mort_aj_cens:pct_desp_recp_saude_mun, 
         desp_tot_saude_pc_mun_def,
         desp_recp_saude_pc_mun_def, 
         ideb_5ano:pop, 
         pct_pop_fem:pct_pop_masc) %>% 
  inner_join(pop_municipio, by = "nomemun") %>% 
  distinct()

Correlação ideb e vacina

# library(esquisse)
# library(ggpmisc)
cor_vac_ideb = dados_vacinas %>% 
  filter(!ano == ymd("2022-01-01")) %>% 
  select(c("nomemun", "ideb_5ano", "ideb_9ano", "cob_vac_bcg", "ano", "porte_municipio", "pop")) %>% 
  ggplot(.) +
  aes(x = ideb_5ano, y = cob_vac_bcg) +
  # geom_point(color = "black", 
  #            shape = "circle") +
  geom_point(aes(fill = porte_municipio, 
                 colour = porte_municipio,
                 size = pop), 
             shape = "circle") +
  geom_smooth(method = "lm",
              alpha = 0.1,
              se = T) +
  scale_y_continuous(expand = expansion(add = 20)) +
      scale_size_continuous(range = c(2,8)) +
  stat_correlation(method = "pearson",
                   size = 3.5,
                   label.x = 5, 
                   label.y = 120,
                   mapping = use_label(c("R", "R2", "P", "n"))) +
  theme(legend.position = "right",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal da BCG e IDEB 5o ano, municípios do estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Municipio",
       y = "Cobertura vacinal (%)",
       x = "IDEB") +
  facet_wrap(~ano, scales = "free", nrow = 2)

cor_vac_ideb

cor_vac_ideb %>% 
  ggsave(file = here("Figuras", "Municipios_2010_2020_BCG_vs_IDEB.png"), 
         width = 20,
         height = 10)

IDHM e vacina

cor_vac_idhm = dados_vacinas %>% 
  filter(!ano == ymd("2022-01-01")) %>% 
  select(c("nomemun", "idhm", "cob_vac_bcg", "ano", "porte_municipio", "pop")) %>% 
  ggplot(.) +
  aes(x = idhm, y = cob_vac_bcg) +
  # geom_point(color = "black", 
  #            shape = "circle") +
  geom_point(aes(fill = porte_municipio, 
                 colour = porte_municipio,
                 size = pop), 
             shape = "circle") +
  geom_smooth(method = "lm",
              alpha = 0.1,
              se = T) +
  scale_y_continuous(expand = expansion(add = 20)) +
      scale_size_continuous(range = c(2,8)) +
  stat_correlation(method = "pearson",
                   size = 3.5,
                   label.x = 5, 
                   label.y = 120,
                   mapping = use_label(c("R", "P"))) +
  theme(legend.position = "right",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal da BCG e IDH municipal, municípios do estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Porte do município",
       y = "Cobertura vacinal (%)",
       x = "IDHM") +
  facet_wrap(~ano, scales = "free", nrow = 2)

cor_vac_idhm

cor_vac_idhm %>% 
  ggsave(file = here("Figuras", "Municípios_2010_2020_BCG_vs_IDHM.png"), 
         width = 20, height = 10)

Taxa de medicos e vacina

variavel_interesse = "tx_med"

cor_vac_txmed = dados_vacinas %>% 
  filter(!ano == ymd("2022-01-01")) %>% 
  select(c("nomemun", variavel_interesse, "cob_vac_bcg", "ano", "porte_municipio", "pop")) %>% 
  ggplot(.) +
  aes(x = tx_med, y = cob_vac_bcg) +
  geom_point(aes(fill = porte_municipio, 
                 colour = porte_municipio,
                 size = pop), 
             shape = "circle") +
  geom_smooth(method = "lm",
              alpha = 0.1,
              se = T) +
  #scale_y_continuous(expand = expansion(add = 10)) +
  scale_size_continuous(range = c(2,8)) +
  stat_correlation(method = "pearson",
                   size = 3.5,
                   label.x = 5, 
                   label.y = 50,
                   mapping = use_label(c("R", "P"))) +
  theme(legend.position = "right",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal da BCG e Taxa de medicos (CENSO), municípios do estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Porte municipal",
       y = "Cobertura vainal",
       x = "Taxa de medicos (1.000 hab)") +
  facet_wrap(~ano, scales = "free", nrow = 2)

cor_vac_txmed

cor_vac_txmed %>% 
  ggsave(file = here("Figuras", "Municípios_2010_2020_BCG_vs_Taxa_Medicos.png"), 
         width = 20, height = 10)

variavel_interesse = "tx_med"

cor_vac_txmed = dados_vacinas %>% 
  filter(!ano == ymd("2022-01-01")) %>% 
  select(c("nomemun", variavel_interesse, "cob_vac_bcg", "ano", "porte_municipio", "pop")) %>% 
  mutate(tx_med = scale(tx_med)) %>% 
  ggplot(.) +
  aes(x = tx_med, y = cob_vac_bcg) +
  geom_point(aes(fill = porte_municipio, 
                 colour = porte_municipio,
                 size = pop), 
             shape = "circle") +
  geom_smooth(method = "lm",
              alpha = 0.1,
              se = T) +
  #scale_y_continuous(expand = expansion(add = 10)) +
  scale_size_continuous(range = c(2,8)) +
  stat_correlation(method = "spearman",
                   size = 3.5,
                   label.x = 5, 
                   label.y = 50,
                   mapping = use_label(c("R", "P"))) +
  theme(legend.position = "right",
        text = element_text(size = 15, color = "black"),
        plot.title = element_text(size = 25),
        plot.subtitle = element_text(size = 20),
        plot.caption = element_text(size = 15),
        plot.margin = unit(c(1,1,1,1), "cm"),
        legend.text = element_text(size = 15),
        strip.text.x = element_text(size = 15, face = "plain"))  +
  labs(title = "Cobertura vacinal da BCG e Taxa normalizada de medicos (CENSO), municípios do estado de São Paulo",
       subtitle = "Período de 2010 a 2022",
       caption = "Fonte: IEPS, PNI, TabNet/DATASUS",
       color = "Porte municipal",
       y = "Cobertura vacinal",
       x = "Taxa de medicos (1.000 hab, normalizado)") +
  facet_wrap(~ano, scales = "free", nrow = 2)

cor_vac_txmed

cor_vac_txmed %>% 
  ggsave(file = here("Figuras", "Municípios_2010_2020_BCG_vs_Taxa_Medicos_normalizada.png"), 
         width = 20, height = 10)

distributions_selec = IEPS_Completo_Regiao_SP %>% 
  filter(ano == ymd("2020-01-01")) %>%
  drop_na(valor) %>% 
  group_by(variavel) %>% 
  mutate(valor = if_else(str_detect(variavel, "pop"), log10(valor), valor),
         valor = if_else(!str_detect(variavel, "pop"), scale(valor), valor)) %>% 
  ggplot() +
  aes(x = valor) +
  geom_density(fill = "#65ADC2",
               color = "#65ADC2") +
  facet_wrap(~variavel, scales = "free")

distributions_selec %>% 
  ggsave(file = here("Figuras", "Região_2020_distribuicao_var_selecionadas.png"), width = 20, height = 15)

QQPLOT

IEPS_Completo_Regiao_SP %>% 
  filter(ano == ymd("2020-01-01")) %>%
  drop_na(valor) %>% 
  group_by(variavel) %>% 
  mutate(valor = if_else(str_detect(variavel, "pop"), log10(valor), valor),
         valor = if_else(!str_detect(variavel, "pop"), scale(valor), valor)) %>% 
  filter(variavel == "cob_vac_bcg") %>% 
  ggplot(mapping = aes(sample = valor)) +
    stat_qq_band() +
    stat_qq_line() +
    stat_qq_point(fill = "black",
                  color = "black") +
    labs(x = "Theoretical Quantiles", y = "Sample Quantiles")